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ABSTRACT 

Spectral Energy Distribution (SED) fitting in the far-infrared (FIR) is greatly limited 
by a dearth of data and an excess of free parameters— from galaxies' dust composi- 
tion, temperature, mass, orientation, opacity, to heating from Active Galactic Nuclei 
(AGN). This paper presents a simple FIR SED fitting technique joining a modified, 
single dust temperature greybody, representing the reprocessed starburst emission in 
the whole galaxy, to a mid-infrared powerlaw, which approximates hot-dust emission 
from AGN heating or clumpy, hot starbursting regions. This FIR SED can be used to 
measure infrared luminosities, dust temperatures and dust masses for both local and 
high- 2: galaxies with 3 to 10-1- FIR photometric measurements. While the fitting tech- 
nique does not model emission from polycyclic aromatic hydrocarbons (PAHs) in the 
mid-infrared, the impact of PAH features on integrated FIR properties is negligible 
when compared to the bulk emission at longer wavelengths. 

This fitting method is compared to infrared template SEDs in the literature using 
photometric data on 65 local luminous and ultraluminous infrared galaxies, (U)LIRGs. 
Despite relying only on 2-4 free parameters, the coupled greybody/powerlaw SED fit- 
ting described here produces better fits to photometric measurements than best-fit 
literature template SEDs (with residuals a factor of ~2 lower). A mean emissivity 
index of /3=1.60±0.38 and mid-infrared powerlaw slope of a=2.0±0.5 is measured; 
the former agrees with the widely presumed emissivity index of (5—1.5 and the latter 
is indicative of an optically-thin dust medium with a shallow radial density profile, 
Adopting characteristic dust temperature as the inverse wavelength where 
the SED peaks, dust temperatures '--^25-45K are measured for local (U)LIRGs, ^--^5- 
15 K colder than previous estimates using only simple greybodies. This comparative 
study highlights the impact of SED fitting assumptions on the measurement of physi- 
cal properties such as infrared luminosity (and thereby infrared-based star formation 
rate), dust temperature and dust mass, for both local and high-redshift galaxies. 

Key words: galaxies: evolution — galaxies: high-redshift — galaxies: infrared — 
galaxies: starbursts 



1 INTRODUCTION 

Modeling galaxies' multiwavelength emission has become 
a sopliisticated effort of extragalactic astronomy. Spectral 
energy distribution (SED) templates, generated by mod- 
eling galaxies' stellar populations and radiation, are used 
prolifically to derive stellar masses, extinction corrections, 
and stellar ages using broadband photometry from the rest- 
frame ultraviolet to infrared wavelengths. They are also 
commonly used to constrain redshifts photometrically (e.g. 
Bolzonella et al.|2000l. The population synthesis-generated 



SEDs used to fit short wavelength data (A < 8fim) are 
complex ([Bruzual fc Charlot"'2003'; Maraston "SooF). They 



depend on the initial mass function (IMF; Salpeter 1955 



Kroupa 2001 Chabrier 20031, metallicity, stellar age, and 
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starbursting timescale — duration, frequency, and strength. 
Although this gives rise to many free parameters in the mod- 
els, the slew of broadband filters in the optical and near- 
infrared make this detailed SED fitting possible, even with 
the effects of dust obscuration/attenuation taken into ac- 
count (e.g. ICalzetti et al.||1994| |Calzetti|[200T| . Follow-up 
spectral observations in the optical and near-IR often con- 
firm good fits to broadband photometry and accurate stellar 
population modeling. 
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The dawn of new infrared observing facilities— from the 
Herschel Space Observatory, the Atacama Large Milhmeter 
Array (ALMA), to the Scuba2 instrument on the James 
Clerk Maxwell Telescope (JCMT)— has triggered a wave of 
interest in extending the use of these template SED libraries 
to the far-infrared (FIR; ~8-1000/im rest-frame), for exam- 
ple, in generating far-infrared photometric redshift estimates 
| |Ros eboom ct al. 2012). Unfortunately, modeling the dust 
emission from galaxies is just as complex in the far-infrared 
as it is in the optical since free parameters include dust 
composition, dust grain type, galaxy structure, orientation, 
active galactic nuclei (AGN) heating, emissivity, and optical 
depth. However, in contrast to the optical and near-infrared, 
far-infrared observations are plagued with a dearth of data. 
Where there might be 10-20 broadbands in the optical/near- 
IR, there are at most ~8 bands in the FIR (most galaxies 
having data in 3 bands or less), all of which suffer from the 
increased beamsize of single-dish FIR observations, thus in- 
creasing uncertainty on measured flux. 



Nevertheless, detailed radiative transfer models and em- 
pirical template libraries have been devised in recent years to 
model the dust infrared emission of stars, molecular clouds 
and starburst galaxies over a wide range of bolometric lumi- 
nosities ( jSilva et a l."1998': ' Chary fc E lbaz"2001': D ale et al 
Helou 20 02, {Kbei 



20011 Pale 



Wandclt 2002j ISieben- 



morgen fc Kriigcl 2007 Draine fc Li 20071. These mod 



els are then often used as a basis on which to measure 
the fundamental parameters of observed starburst galaxies 
with constrained photometric measurements, both at low 
2009] [AFmus et al.|[2009l [Chapin 



redshift (e.gJRieke et al. 



et al. 



et al. 



2009 



2002 



U et al 



2012 1 and at high redshift (e.g. [Blain 



Chapman etaL][2004l |Pope et al.|[2008l [Swm 



bank et al.| 2010 1. Particularly for galaxies with high star 
formation rates, the 'UV chimney' argument (Neufe ld|1991 1 
can be used to relate the modeled dust and molecular cloud 
structure back to the Ly-a escape fraction; Ly-a suffers 
from less attenuation than UV continuum photons despite 
very large reservoirs of dust and gas. 



This paper investigates the use of these far-infrared 
template SEDs to determine fundamental galaxy properties 
such as infrared luminosity. Lib., characteristic dust tem- 
perature, Tdust, and dust mass, Mdust- Section |2] presents a 
simple method for representing galaxies' FIR emission as a 
coupled modified greybody plus a mid-infrared (MIR) pow- 
erlaw, which can constrain these fundamental IR-derived 
properties quite well for a wide range of galaxies. Sec- 
tion [3| compare s these fits to the SED template fits of 



Char>rfc Elbaz| ([200l|), [Dale fc Helou| p002t and [Sieben- 



morgen & Kriigel (20071, using the photometry of lo- 



cal luminous infrared galaxies and ultraluminous infrared 
galaxies (LIRGs/ULIRGs) from |U et al.| (|20T2|. Section 
|4] discusses derived quantities of FIR SED fits, and sec- 
tion [5] concludes. Throughout, a A CDM cosmol ogy with 
-ffo=71 kms"'^ Mpc~^ and r2m=0.27 is assumed (Hinshaw 
et al.|2009| ). 



2 SED FITTING TECHNIQUES 

2.1 Coupled Greybody/Powerlaw Fitting 

2.1.1 Method 

Deriving the fundamental physical properties of IR- 
luminous galaxies can be as simple as assuming an isotropi- 
cally emitting blackbody. This is represented as the Planck 
function, Bv{T) (e.g. in units of ergs~^ cm~^ A~^), and is 
only dependent on dust temperature T. However, if the vari- 
ation in opacity (e.g. assuming a screen of dust without 
scattering) and source emissivity is accounted for (the fact 
that very few sources are perfectly non-refiective) , the flux 
density at rest-frame frequency v is then represented by a 
modified black body (i.e. "greybody") of the form 



S{iy) oc (1 - 

where S{h' 
depth is t{i' 



(1 



(1) 



is in units of ergs cm Hz 
and fitted as t{u) — [v/uoY 



or Jy. Optical 
, where is the 



frequency where optical depth equals unity ( Draine 2006 1 



and /3 represents emissivity, or the spectral emissivity index. 



See Kovacs et al. ( 2010 1 for a thorough discussion of the im- 



pact on /3. The value of P is largely assumed to be 1.5 (and 



usually ranges 1-2, Hildebrand 1983 1, although this could be 



a result of the original wavelengths for which data were gath- 
ered on local starbursting samples ( Dunne fc Eales^^200l| ). 
Some recent work points to a wider range of /3 values be- 
tween 1-2.5 (e.g. [Chapin et al.|201l] [Casey et al.|2011[ ). The 
theoretically expected value of vq is 3THz (i.e. Ao=100/im), 
although this value is unconstrained by data (see discussion 
Conley et al.[[20TT l. In the optically thin case, the term 



•') reduces to v'^ , and the flux density simplifies to 



(2) 



The normal range of dust temperatures expected for a 
galaxy's interstellar medium (ISM) heated only by star for- 
mation ranges ~20-60 K. When fitted to a greybody (as in 
Equations [1] or [2|, most galaxies have a notable flux density 
excess at wavelengths shortward of «50/im (see Figure [l] 
panels A and B). This mid-infrared excess is due to a com- 
bination of hotter-dust subcomponents (where dust is more 
compact) or dust heated by an AGN, and an optically thin 
medium by which the higher frequency radiation can escape. 
The disconnect between observed mid-infrared luminosities 
and predicted Wein-tail luminosities has been studied for 
quite some time; a quite thorough discussion of dust clouds' 
opacity, radial density distributions, and dust mass coeffi- 
cients (k,v), impact on observed SED is given in [Scoville[ 



fc Kwan (19761; note in particular the difference between 
optically thin and optically thick models in the near- and 
mid-infrared regime. 

Non-thermal emission from polycyclic aromatic hydro- 
carbons (PAH) in starburst galaxies or silicate absorption 
at 9.7^m can contribute to the SED shape in this regime as 
well; however, their net effect on integrated IR luminosity, 
differing from a simple mid- infrared powerlaw, is <10% in 
most cases (See 



2.21 



While the hot-dust component might be made up of 
several subcomponents of different warm temperatures, the 
cold dust modified greybody still dominates the bulk of 




Figure 1. Illustration of FIR SED fitting techniques using the example local galaxy UGC02369. Panel (A) highlights the four IRAS-hand 
photometric points in black (the more recent FIR photometric measurements in gray). The original computed FIR luminosity for this 
source was computed via Eqjs] The original dust temperature for the source (Perault 1987) was found by fitting a single temperature 
greybody with fixed /3=1.5 to the IRAS bands (magenta; the same best-fit SED assuming optically thin conditions is shown in dashed 
magenta). Panel (B) shows a slightly improved SED fit, with a single temperature greybody fit to all of the FIR photometric bands. 
Emissivity, /3, is added as a free parameter (black) and measured as /3 =0.5, which is similar to adding a very cold greybody component 
to the warm component dominating the emission, both with fixed /3 =1.5 (e.g. 60 K and 15 K SED shown). Both the fits in panels (A) 
and (B) have mid-infrared excesses and are poor fits to a single temperature greybody. Panel (C) shows an improved fit, using four 
temperature greybodies with fixed /3=1.5. The fit to the data is much better than (A) or (B), although the number of free parameters 
shoots up to 9. Panel (D) illustrates the SED fitting technique described in § |2.1[ a composite mid-infrared powerlaw plus cold-dust single 
temperature greybody. The mid-infrared powerlaw is a good approximation for the composite of warm dust giving rise to the mid-IR 
excess. All fits here assume general opacity conditions; in the optically thin case, the fits alter slightly, but produce the same Lir and 
Tdust to within '~1%. 



the total infrared emission when integrated. The net sum 
is an SED which can be approximated as a powerlaw 
in the mid-infrared, with intensity dropping with decreas- 
ing wavelength, and a single temperature greybody fit in 
the far-infrared. The greybody dominates at wavelengths 
>50/im whereas the mid-infrared powerlaw dominates at 
wavelengths <50/im. This can be analytically approximated 



20111, P represents the emissivity, a represents the slope of 



the mid-infrared powerlaw component, and Ac is the wave- 
length where the mid-infrared powerlaw turns over and no 
longer dominates the emission. This simplifies to 



•Sot (A) = A'"bb,ot - 



Jli +iV,A"e-(^' 



(4) 



^(A) = iVbb 



(1 



^hc/>kT . 



+ iVpi A'^e-f^' 



(3) 



where S(A) is in units of Jy, T is the galaxy's characteristic 
"cold" dust temperature (in other words, the dust temper- 
ature dominating most of the infrared luminosity and dust 
mass), Ao is the wavelength at which optical depth is unity 



(taken here to be fixed at Ao=200^m as in Conley et al 



in the optically thin case. Note that several works have rec- 
ognized the utility of adding a mid-infrared powerlaw com- 



ponents to simple greybody fits (e.g. Younger et al. 2009 1 



but tend to fit the two components separately: first the grey- 
body followed by the mid-infrared powerlaw. Coupling the 
two together and fitting simultaneously makes it possible for 
a more accurate fit to be made to systems with fewer FIR 
photometric datapoints. 
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Figure 2. At left, a comparison of tlie T temperature parameter of tlie fitted greybody against tiie "peak wavelength" Tj, which is 
proportional to the inverse peak in Si^. This varies by fitting method mostly due to slight variations in adopted opacity, which shifts 
the SED peak slightly while keeping Ljpj^ roughly constant. At right, a comparison Tj, the "peak wavelength" temperature, to IRAS 
[60]-[100] color, which has often been used as a proxy for dust temperature. Over-plotted are various model template values, from |Chary| 
|fc Elbaz| pool} , [Dale fc Helou| pOOl] !, [Siebenmorgen fc Kriigel| | [2507l l and [Draine "fcXIl p007| l . 



Figure[T]illustrates this SED fitting technique relative to 
single temperature greybody fits (as in Eq[T]) for an example 
local LIRG, UGC 02369. This highlights the inconsistency 
of single temperature greybodies with FIR photometry of 
real galaxies and the usefulness of an SED fitting technique 
which simultaneously fits the cold-dust, long- wavelength 
greybody and mid- infrared excess powerlaw component. 

Another common way to consider the 'warm dust' com- 
ponents is as a continuous powerlaw distribution of tem- 
perature components, e.g. dMd/dT oc T"^, as described in 
Kovacs et al. (20101. The index j effectively represents a; 



the values 7=7.2 and 7=6.7 (both assumed in the Kovacs 
et al.|[2010 analysis for different samples) correspond to a 
values of ~2.9 and 2.6 respectively. Kovacs et al. (20101 
points out that 7~6.5-7.5 is expected for sources in a dif- 
fuse, star forming medium, and 7~4-5 for a dense medium 
from the heating source (e.g. a dusty torus around an AGN). 
The fiux density function is then given as 5y(T^_>cx>) (r) ~ 
(7-l)rT-^ Su,T{T)T-'<dT. This integral must be solved 
numerically, but leads to an SED of shape similar to those 
produced by Equations [3] and [4] particularly Equation [3] 
The differences between the temperature powerlaw integral 
fit and the analytic approximation in Eq [3] is that the lat- 
ter provides a bit more flexibility in that the relationship 
between normalization factors can be adjusted. The former 
has a clear, robust motivation, but is also a computation- 
ally expensive algorithm. As photometry for infrared sources 
improve and spectroscopy becomes more widely available, 
a minority of galaxies will likely have unusual SED shapes 
not well described by dM^/dT oc T""*", in which case a more 
general form of Eq[3] might be appropriate. 

Note that another popular infrared SED fitting method 
joins a greybody to a mid-infrared powerlaw using a piece- 



straightforward way of generating SEDs of the desired shape 
(panel D of Figure[T]). In many instances, this method might 
be preferred over this paper's fitting techniques (e.g. when 
generating hypothetical SEDs for testing selection, or gener- 
ating fits for galaxies with only one FIR photometric point). 
However, one disadvantage of this method is its piece-wise 
nature, making it difficult to fit data across the whole in- 
frared range simultaneously and quantify the errors on each 
parameter. Having an analytic approximation for this func- 
tional form (e.g. Equations [s] and |4| makes error propaga- 
tion and multiple-datapoint fitting much more straightfor- 
ward. 



2.1.2 Adjusting the 2-4 free parameters 

There are six total parameters in this fit: the greybody nor- 
malization (A'^bb), the powerlaw normalization (Afpi), grey- 
body temperature (T), emissivity index (/3), mid-infrared 
powerlaw slope (a), and mid-infrared turn-over wavelength 
(Ac). Since the turnover wavelength of the mid-infrared pow- 
erlaw would depend on the turnover point of the cold dust 
greybody component and A^pi determines the flux scaling of 
the powerlaw term relative to the greybody, neither param- 
eter should not be considered 'free' in the sense of the other 
four. Both Ac and TVpi are tied to the best-fit values of A'bb, 
T and a such that the total SED resembles both the lKovacsl 



et al. (20101 powerlaw temperature distribution SED and 



wise technique (Blain et al. 2002 20031, where the tran- 



sition point between greybody and powerlaw is defined by 
d/d\ {S)=a, i.e. the gradients are equal. This is the most 



the Blain et al. (20031 piece- wise matched-gradient SED. 

The turnover wavelength, Ac is set to 3/4 the wave- 
length where the gradient of the greybody is a, as in 
Blain etal.. Ac is both a function of a, T and the adopted 
opacity model, and can be approximated as Ac = (T x (ai -f 
a2a))~^ (optically thin) and Ac = ((bi + 620)^^ + (bg + 
hio) X r)~^ (general opacity). The values of the coefficients 
are given in Table [l] Note that the factor of 3/4 is incorpo- 
rated so that the juncture of the powerlaw and greybody is 
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most smooth (i.e. it does not have a physical interpretation, 
although is related to the falloff rate of ) . The coefRcient 
of the powerlaw term, A^pi, is tied to the normalization of 
the greybody term and solves the condition A^piA" =S^(Ac), 
where S^, is given by Equations [l] or [2] Having fixed these 
two parameters reduces the number of free parameters to 
2-4. 

Depending on the amount of FIR photometric data 
available, further constraints can be made to reduce the 
number of free parameters in the fit described in EquationjS] 
or Equation |4] from four (A'^bb, T, /3, and a) to two (A^'bb 
and T). The number of free parameters should never ex- 
ceed the number of independent data points minus one. The 
emissivity jS varies from 1-2.5 in the literature for individ- 
ual sources, although the vast majority of works assume a 
fixed value of;3=1 .5 ( [Chapman et al.|2005|rPope et al.|2"006 



Casey et al.]|2009 among others). Varying j3 does not have 



a very significant impact on FIR luminosity or dust tem- 
perature, but does have a strong impact on the slope of the 
Rayleigh- Jeans tail at rest-frame A >200//m. If there are 
>3 independent photometric points at A >200/im, then /3 
can be constrained with the fit; otherwise, the fixed value of 
/3=1.5 is suggested. The mid-infrared powerlaw slope can be 
constrained similarly; if >3 photometric points are available 
at rest-frame A <70/im, then a can be measured. Other- 
wise, a fixed value of a=2.Q is consistent with most sources 
and is directly comparable to the mid-infrared portion of 
SED templates presented in the next section. Typical values 
of a range from 0.5 (nearly flat, consistent with quite a bit 
of warm dust) to 5.5 (very steep, consistent with very little 
warm dust). Note that values of a <1 lead to a divergent 
luminosity in the near-infrared, so some additional cutoff 
should be placed at short wavelengths to avoid a nonphysi- 
cal interpretation. 



2.1.3 Defining Lib. and Tdust 

The normalization of the fits is governed by the infrared lu- 
minosity, LiR, whose integration limits have varied through- 
out the literature. Lir is intended to represent the bulk of a 
galaxy's dust emission. For galaxies which are very infrared- 
bright (LiR >10^^ L©), Lm_ is a proxy for bolometric lumi 
nosity. Lir is most often taken from 8-1000 /im (e.g. Kenni 



cutt 



1998 1. Lir also goes by different names in the literature: 
Ltir ("total IR") or Lfir ("far IR") and can be integrated 
from 3-1100 /xm (e.g . [Chapman et al.|2003[), 40-120 /^m (e.g . 
Younger et al.[2009| ), or 40-1000 ^m (e.g. [Conley et al.[2011| ). 
The rest of the paper uses the standard 8-1000/.im integra- 
tion limits for the sake of consistency with most of the liter- 
ature; however, if only the cold-dust, star-formation domi- 
nated infrared luminosity is desired, more appropriate limits 
would be 40-1000/im. 

In the initial years of characterizing /7?j4S'-detected IR- 
luminous galaxies, the integrated IR fiux was approximated 
by: 

-FpiR(40-500) ~ 1.26 X 10"" {2.58/60 + /100} (5) 

in Wm"'^ (and where flux densities, /, are given in Jy, see 
the review in [Sanders fc Mirabel[[1996[ ). Including all IRAS 
bands and integrating over a wider wavelength range, 8- 
lOOO/xm, this becomes: 



-FlR(8-1000) — 1.8x 10-'*{13.48 /12+5.I6 /25+2.58 /eo+Jioo} 

Luminosities are then LiR=47r D|,Fir. Section [3] discusses 
the accuracy of this widely accepted luminosity approxima- 
tion to the |U et al.| sample. 

Note also that the dust temperature T as given in Equa- 
tions [l||4| represents the intrinsic galaxy temperature, which 
is different than the inverse "peak wavelength" temperature 
as measured by Wein's displacement law, i.e. \max=b/Td 
(where fo=2.898xl0''/imK), which only applies to perfect 
blackbodies. Figure[2|shows the intrinsic dust temperature T 
against peak wavelength dust temperature, Td for the above 
formulations. Figure [2] also shows Td against IRAS color, 
log{Sfio/Sioo), for Equations [lj|4] and the IR template SEDs 
discussed in the next section. 

Far-infrared color is often taken as a proxy for dust 
temperature. Works in the literature vary the use of differ- 
ent fltting methods, some using IR template fits, some using 
optically thin greybodies, some including a warm dust com- 
ponent and some not. The convention of referring to a single 
dust temperature is meaningless unless it refers to the same 
measurement made with consistent methodology. However, 
the Wein's displacement law value of dust temperature, Td, 
does not vary greatly between fitting techniques, so it is in- 
sensitive to biases in fitting methods, making it a good diag- 
nostic of the characteristic dust temperature for the whole 
system. Throughout the rest of the paper, Td is used to 
represent galaxies' dust temperatures and this method also 
permits fair comparisons across literature sources. 



2.2 IR Template SEDs 

Template spectral energy distributions are taken from four 
sources, all of which are widely used in the literature to de- 
rive values such as Lfir (thus star formation rate, SFR), 
Tdust, and Mduat for nearby and high- 2; galaxies. The first 



105 template SEDs are from Chary & Elbaz (20011, whose 
models are developed to represent existing data from ~0.4- 
850/im on nearby galaxies (and which vary as a function 
of LpiR, Tdust and PAH strength). We also use the 64 phe- 



nomenological models of Dale et al. ( 2001 1, supplemented by 

The SEDs 



FIR/submillimeter data in Dale & Helou 



(20021 



from |Chary fc Elbaz[ ( |200T| ) aiid fDale fc Heiou| ( |2002|| are the 



models most frequently used in the literature (e.g. Chapman 



et al.|2005|[Chapin et al".t2009, and many others). However, 
the radiative transfer model SEDs from [Siebenmorgen fc| 
Kriigel (20071 are also included, and their parameter space 



spans much more diverse SED types (by varying nuclear 
radius, visual extinction, ratio of OB stellar luminosity to 



total, and dust density). While the Siebenmorgen & Kriigel 
( |2007l ) work contains ~7000 SEDs, we limit ourselves to the 
120 most luminous SEDs, consistent with ULIRG infrared 
luminosities. Finally, we include 25 of the most IR-luminous 
models from the PAH-motivated work of Draine & Li ( 2007| 
which have been used for high- 2; samples (e.g. see ,Marsden| 
et al.|2011[ |. 

The primary observable difference between the in- 
frared template SEDs in this section and the coupled grey- 
body/powerlaw is the inclusion of the PAH emission fea- 
tures and silicate absorption features. Readers interested in 
the relation between PAH emission and FIR dust emission 
should clearly use template SEDs instead of simple grey- 



6 C. M. Casey 



5(A) = ATbb 



Free Parameters: 



Nbb, T, (13, a) 



Emissivity 

Mid-IR Powerlaw Slope 
Wavelength where opacity is unity 
Powerlaw turnover wavelength 



13 1.60±0.38 
a 2.0±0.5 
Ao =200 Atm 

Ac =3/4 L(a,T) 
where L(a,T) = ((61+620)"^ + {bs+bi a)x T)-'^ 

bi 26.68 

62 6.246 

63 1.905x10-* 

64 7.243xl0-'5 
Normalization of powerlaw term A^pi 



(l-e-(^0/^c)'^);^-3 

Characteristic Dust Temperature T) 1/ Apeak 



Table 1. Adopted Best-Fit SED fitting equation and measured 
parameters based on the GOALS sample. The measurement of 
/3 is based on 48 out of the 65 galaxies which have 850— 1.1mm 
data, while the measurement of a is based on the full sample of 
65 (U)LIRGs. 



body/powerlaw fits. However, in terms of the basic phys- 
ical properties extracted from FIR SEDs— Lir, Tdust, and 
MduBt , there is not a significant detriment to using the grey- 
body /powerlaw fitting technique. The PAH emission lines at 
7.7/im, 11.2/xm, and 12.8/im contribute to the integrated IR 
luminosity on the order of 5%, although their contribution 
is often negated by the presence of the 9.7fim absorption 
feature. 

The mid-infrared spectra of dusty starbursts have been 
shown to vary substantially at both \ow-z and high-z 
I Brandl et al]|2006| [Pope et al!]|2008| [Menendez-Delmestre 
|et al.||2009[ ). Without a direct mid-infrared spectrum, it is 
impossible to know whether or not PAHs or absorption fea- 
tures are contributing significantly to broadband photomet- 
ric measurements. At certain rest-frame wavelengths, FIR 
photometric measurements are likely to deviate positively 
or negatively from a mid-infrared powerlaw (e.g. at ~8/im 
in the former case and lOfim in the latter case). In those 
cases, infrared template SEDs should be used to more accu- 
rately determine the structure of mid-infrared emission. 



3 SED COMPARISONS WITH DATA 

Testing and comparing the different IR Template SEDs from 
Section |2.2| and SED fitting methods from Section |2.1| re- 
quires a sample of well-constrained FIR-bright galaxies with 
accurate FIR photometric measurements. To date, the most 
extensively imaged IR galaxies— besides the handful of high- 



z brightly lensed sources (e.g. the "cosmic eyelash," Swin- 
bank et al.|20T0 | —are from the local /7?^5'-selected Revised 
Bright Galaxy Sample (RBGS; Sanders et al.|2003 l, partic- 
ularly the subset which form part of the Great-Origins All 
Sky LIRG Survey (GOALS ; |A"rmus et al.|2069j ). The revised 
standard-aperture photometry for 65 galaxies in the GOALS 



sample is summarized and presented in U et al. (20121, and 



used herein at rest-frame wavelengths 5-2000/im. 

Modified greybodies plus mid-infrared powerlaw SEDs 



are fitted to data as described in Section [2. 1[ specifically to 
Equation [3] taking both flux density uncertainties and non- 
detection upper limits into accounlQ The general-opacity 
fits are chosen for this local sample over the optically-thin 
fits due to a subtle difference in SED shape around rest- 
frame wavelengths ~20-40/im. 

For sources with more than three long-wavelength data 
points, P is kept as a free parameter; otherwise it is fixed 
to 13=1.5. Mid-infrared powerlaw slope a is allowed to vary 
since the number of short-wavelength points (<70^m) is al- 
ways more than three. The range of a values is found to 
vary between a=0.5-5.5. Dust temperature and luminosity 
are left to vary. Template IR SEDs from Section |2.2| are 



normalization (i.e. Nbb ocLir) is a free parameter. 

Figure |3] shows several randomly-selected examples of 
best-fit SEDs for GOALS sources. The various best-fit IR 
templates are colored, while the fitted SED is shown in 
black. The underlying greybody function to the fitted SED 
is shown as a dotted line. Since >8 FIR photometric points 
exist for the GOALS sample, a dummy SED can also be 
constructed by linearly extrapolating between photometric 
points (shown in gray). While this is a rough approximation, 
the advantage of this extrapolation is that this dummy SED 
makes no intrinsic assumptions as to SED shape. The sub- 
set of sources shown in Figure [3] represents a wide variety of 
templates and fitted SEDs and is representative of the fits 
to the whole [U et al. ( 2012[ ) sample. 

The residuals of the fits with respect to data are shown 
in Figure [4] The residual is computed as the difference in 
logS*^ between the data points and the best fit SED at each 
of the following wavelengths: 12/^m, 25/im, 60/im, lOO/im, 
and 850/im. Figure [4] shows the distribution in the residuals 
at each wavelength for the |U et al.| GOALS sample for each 
fit. Black represents the residuals from the fitting method 
described in Section |2.1| and the residuals of the IR tem- 
plates are colored as in Figure [S] At each wavelength, the 
mean and standard deviation of the residual distribution is 
given in the top left. 

At all wavelengths, the fitted SED from Equation [3] 
is statistically a better fit to the data than any IR tem- 
plates, despite having fewer free parameters than the SED 
templates; a summary of derived quantities from this fit is 
given in Table[l] However, some templates are more accurate 
than others; the [Siebenmorgen fc Kriigel| templates provide 
the best fit of the templates, significantly better than both 
|Chary fc Elbaz| and [Dale fc Helou| templates. However, all 
seem to suffer at short wavelengths where there is natu- 
rally no flexibility built into the models. This demonstrates 
that galaxies' mid-infrared properties are not tightly cor- 
related with their far-infrared properties, as perhaps previ- 



ously thought (since 24/im is often used to infer Lir, e.g. Le 



Floc'h et al.|2005 on Spitzer samples, among many others) 



Another reason the [Siebenmorgen fc Kriigel| templates 
have lower residuals than the other two template libraries 
is that they contain more templates with which to compare 
the data. Even at the highest luminosity end (Lir >10^^Lq), 



^ An IDL function cmcirsed.pro, which can be used to fit an 
SED of the form in Eq|3]or Eq^to real data, is publicly available 
at www. if a.hawaii . edu/^cmcasey /research. html. 
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Figure 3. Example best-fit SEDs to a subset of the GOALS sample, whose photometry is presented iii |U et al.| l |2012| |. Photometric 
data points are shown as red diamonds with uncertainties illustrated (sometimes smaller than the size of the data point). The best-fit 
SED fit from Equation |3] is shown in black, with the underlying greybody distribution shown as a dotted black line. The best-fit [Chary| 
I&: Elbaz| template SED is shown in magenta, |Dale fc Helou] SED shown in green, and |Sicbcnmorgcn Sz Kriigel| SED shown in orange. A 
"dummy SED" is also constructed which linearly interpolates between data points in logA-logS^ space (dashed gray line). The 8-lOOO^m 
and 40— 1000/im integration limits are marked as light gray vertical lines for reference. 



there are still 120 template models i n the |Siebenmorgen fc| 
Kriigel ( 2007 1 library. In contrast, the Chary & Elbaz library 
contains 105 template SEDs and [Dale fc Helou| contains 64 
template SEDs. 

This direct comparison between SED fitting and SED 
templates demonstrates that using IR templates will always 
provide a less accurate fit than simple greybody/powerlaw 
fitting. This is not meant to suggest that IR template SEDs 
are not of great use; on the contrary, they provide very good 
constraints on the relationships between FIR emission and 
PAH emission, and sources' dust composition, which simple 
greybody fitting cannot do. However, if the purpose of FIR 
SED fitting is to measure a source's FIR luminosity, dust 
temperature and dust mass accurately, this work demon- 
strates that direct fitting, as in Section [2. 1[ is best. 

3.1 Interpretation of best-fit parameters 

The |U et al.| ( |2012[ ) sample has been used to constrain the 
mean SED parameters j3, the spectral emissivity index, and 
a, the mid-infrared powerlaw slope. How do we interpret the 
measurements /3=1.60±0.38 and a=2.0±0.5? Here I draw on 



the detailed discussions in Scoville & Kwan ( 1976 1 



The measured value of the emissivity index, /?, is spot on 
its presumed value of ~1.5. Although not new information, 
this suggests that even the cold dust component (dominating 
at A>50/im) has a temperature gradient whereby the dust 
farthest from the luminous starburst is slightly colder, which 
deferentially boosts the fiux density at the longest wave- 
lengths. This boosting (above a simple blackbody) becomes 
enhanced when the density distribution falls off slowly, and 
with a steep density drop off, /3 should increase, and the 
far-infrared "slope" should increase. 

The mid-infrared powerlaw slope, a, probes the warmer, 
more-compact dust, and can be used to estimate the radial 
density profile. Our measurement of a=2.0±0.5 agrees with 
earlier findings of a =1.7-2.2 from analysis of IRAS and 
distant galaxies in |Blain et al.| ( |2002[ ). Assuming there are 
sufficient hot dust grains near the interior and the dust is 
optically thin, S\ oc A° where a=2.0±0.5 translates to a 
density profile of y-0 5±0-2^ consistent with r^^^^ . In other 
words the dust density profile is relatively flat and perhaps 
indicative of the diffuse nature of the dust around these star- 
bursts, their distributions perhaps originating from the vio- 
lent interacting nature within the galaxies. 
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Figure 4. The residuals of the fits to the GOALS sample be- 
tween template libraries and the coupled greybody / powerlaw fits 
for all 65 galaxies, for the greybody /powerlaw fits (black), Chary 
Sz Elbaz fits (magenta). Dale & Helou fits (green) and Sieben- 
morgen &; Kriigel fits (orange). The mean and standard devi- 
ations of the residual distributions for each fitting method are 
given in the upper left-hand corner. The residuals are lowest for 
greybody/powerlaw fits, followed by the Siebenmorgen & Kriigel 
fits. 



4 DERIVED QUANTITIES 



Section [s] u sed the GOALS sample ( [Arrnus et al]|2009| |U| 
|et al.|2012[ ) to demonstrate that direct greybody/powerlaw 
SED fitting, as described in Section [2. 1[ is a more accurate 
fit to data than best-fit IR SED template hbraries from the 
hterature. In this section, the derived quantities Lir (lu- 
minosity), Td (dust temperature) and Md (dust mass) are 
compared to see what can be expected using difi'erent SED 
fitting methods. 



4.1 IR Luminosity and Star Formation Rates 

The IR luminosity is most often integrated in the range 
8-lOOOAtm, i.e. Lm = 47rDi J^^°^ S^dv. The range 8- 
1000/jm encompasses some of the most prominent PAH 
emission features in addition to the greybody emission peak. 
The boundaries are stih somewhat arbitrarily drawn, since 
8/im sits awkwardly on top of the 7.7/im PAH emission fea- 




11.5 12.0 
L,R [Lo] IRAS 



Figure 5. The fitted RBGS IR luminosities, "IRAS luminosi- 
ties," are measured as a linear combination of IRAS 12-100/^m 
fiux densities, see Eg [e] or [Sanders Mirabel) l |1996| |. Here IRAS 
luminosities arc plotted against fitted luminosities from the direct 
fitting method (Section |2.1| , and template SED libraries. The 
full RBGS sample of 625 sources from [Sanders et al.| l [2003[ l is 
included, and luminosities refitted using IRAS photometry alone 
and Eq[3| fixing /3=1.5 (scatter is larger for this sample due to 
poorer photometric constraints). At bottom are the residuals, lu- 
minosity against the difference in logLjR. All template libraries 
and fitted SEDs demonstrate that IRAS luminosities are overes- 
timated by ~0.04dex, or ~10%. 



ture, and lOOO/xm splits the cold-dust greybody neither at 
its turnover point, where radio synchrotron emission begins 
to dominate, nor near the peak. On average, the cold-dust 
greybody component— thought to be exclusively heated by 
star-formation processes— comprises 74±11% of iiR(8-iooo)i 
while warm AGN-heated dust and PAH emission compo- 
nents comprise 26±11%. 

Figure [5| shows the 8-1000/im IR luminosities of the 
[U et alT] sample and luminosity residuals when contrast- 
ing methods. The x-axis is the IR luminosity as inferred by 
Equ ation[6[ also as summarized in'Sanders fc Mirabel[ ( [1996[ ) 
and Armus et al. (2009). The whole RBGS sample is also 



included, refitting all luminosities using the method from 
Section [2.1| (due to limited long-wavelength photometry, 
13—1.5 is fixed). Similarly, RBGS luminosities are compared 
to template SED-inferred luminosities. Both the fitted SEDs 
and template SEDs have lower IR luminosities than pre- 
dicted using the RBGS Equ ation (Eq[6| by ~0.02-0.04dex. 
Both [Chary fc Elbaz[ and [Dale fc Helou[ templates are 
0.02±0.07dex lower in luminosity than IRAS predictions, 
while [siebenmorgen fc Kriigel[ templates are 0.03±0.05dex 
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as in Conley et al. (2011 1. Longward of 40^m, the cold-dust 
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Tdust [K] Perault Greybody 

Figure 6. A comparison of dust temperatures from Per- 
ault (1987), described also in [Sanders fc Mirabel] l |l996| , against 
dust temperatures found using methods in SectiorL |2.1| and from 
template SEDs. Colors and symbols are as in Figure [5] In all 
cases, a single greybody fit of the form described in Eq [l] or [2] 
(which was the method used in Perault 1987) overestimates dust 
temperature by 10 K on average. Caution is necessary when using 
dust temperatures from template SEDs since they are quantized. 



lower, and fitted SEDs are 0.04±0.04dex lower. These dif- 
ferences translate to luminosity overpredictions from RBGS 
of 10%, thus overestimated SFRs by 10%. 



Note that the Kennicutt ( 1998 1 scaling relation between 



^iR(8-iooo) a-nd SFR take both the contribution of "infrared 
circus" cold-dust emission and "warm" dust around young 
star forming regions into account and is generated using stel- 



lar synthesis models ( Leitherer & Heckman 1995 1 for con- 



tinuous bursts aged 10-100 Myr. Since the Kennicutt scaling 
is not determined as an empirical relation between, e.g. Ha 
SFR and Lir, there is no need to recalibrate it to correct 
for the overestimated Lir; however, the FIR star formation 
rates of the local (U)LIRG sample should be revised to re- 
flect the systematic offset. 

Regardless of the systematic offset, infrared-based star 
formation rates are only accurate to ~10% mostly due to the 
variation in mid-infrared properties, likely not directly scal- 
ing to the galaxies' star formation rates, since complex pro- 
cesses dominate the mid-infrared. Similarly, slight changes 
in any SED which impacts the mid-infrared (e.g. changes in 
opacity assumptions, assumed a, PAH contributions) will 
impact the derived star formation rate. For this reason it 
is worth emphasizing that the systematic offset is of minor 
significance. Also worth pointing out is that the range 8- 
lOOO/xm is not ideal to generate precise star formation rates. 

Ideally, one would extract the integrated luminosity 
from the cold-dust greybody component over all wavelengths 
where it dominates. This can be estimated precisely if the 
fitting method in Section [2?T] is used, as the greybody compo- 
nent can then be constructed independently after fitting to 
a joint greybody and powerlaw. However, to make fair com- 
parisons with different SED fitting methods, the optimal al- 
ternative limits of integration should be taken as 40-lOOO/im 



greybody dominates at all dust temperatures <100K. 



4.2 Dust Temperature 

Comparing dust temperatures between models first requires 
an understanding of the dust temperature convention for the 



original work on the U et al. |2012) GOALS sample objects. 
Dust temperatures were originally fitted in Perault (1987; 
PhD Thesis) by fitting a single temperature dust emissivity 
model (e oc v~^) to the fiux in all four IRAS bands (briefly 
described in [Sanders fc Mirabel|1996 1. In other words, these 
dust temperatures are fitted to data ranging 12-lQO/im to 
either Eq [l] or Eq [2] Assumptions made with regard to 
opacity are not stated; however, the impact on peak dust 
temperature is minimal. It could be that the dust tempera- 
ture reported in Perault (1987) is indeed T taken from e.g. 
Eq [2] rather than "peak wavelength" dust temperature Td, 
although it is unclear. Regardless, refitting greybodies of ei- 
ther form specified in Eq[l]or [2]to only the four IRAS da,t- 
apoints causes the dust temperatures Td to drop by ~10K. 
This is shown in Figure [6] (and was illustrated clearly in Fig- 
ure [l] panels A and B). The difference between dust tem- 
peratures of simple greybodies and greybodies which incor- 
porate a mid-infrared component has been known (see also 
the recent detailed discussion in Hayward et al. 20121, yet 



few works are explicitly clear on how their dust temperature 
measurements should be interpreted. 

The Perault/RBGS dust temperatures are also found 
to be significantly hotter than dust temperatures measured 
from the SED template libraries. However, any dust tem- 
peratures extracted from the template SEDs should be used 
with great caution, as the dust temperatures are quantized 
in all cases (as was seen at right in Figure [2|. The Dale 
[fc Helou] (12002') templates exhibit the poorest constraints 
on dust temperature, as the templates only have two dust 
temperatures: Td Ri31K and Td Ri51K. The [Chary fc Elbaz[ 
( |2001[ ) templates have four different dust temperatures, at 
Ri23K, 32 K, 41 K and 50 K. The [Siebenmorgen fc K"rugel| 
(2007) templates are quantized on yet finer scales, at 5K in- 
crements between ~24 K and 85 K; however, the fitted SEDs, 
as an analytic fit, are not quantized in Td. For this reason, 
any dust temperature-based results should be based on fits 
of the type described in Section [2. 1[ not on template SEDs 
in the literature. 



4.3 Dust Mass 

Dust mass is related to IR fiux density and dust temperature 
via 



(7) 



where S^ is the fiux density at frequency u, is the dust 
mass absorption coefficient at v, B^(r) is the Planck func- 
tion at temperature T, Md is the total dust mass, and Dl 
is the luminosity distance. The fact that our sources are not 
perfect black bodies is accounted for by the dust mass coef- 
ficient Ki, so that the greybody is effectively represented by 
the product K^BviT), and the sources' luminosity at fre- 
quency V scales as S„/Bi,{T)cc . The dust mass is then 
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Figure 7. Dust mass against S-lOOO/^m IR luminosity and star 
formation rate (SFR) measured from the infrared. Isotlierms are 
marked with dashed lines increasing by 10 K. This emphasizes 
that (a) dust mass is incredibly sensitive to dust temperature 
and that dust mass should only be computed when dust temper- 
ature is adequately constrained, and (b) the lack of correlation 
between integrated IR luminosity and the Rayleigh-Jeans cold- 
dust component of IR SEDs (the best probe of dust mass). The 
latter point points out that using direct scaling relations between, 
e.g. 5g5o and SFR, is not recommended. 



D 



(8) 



As Draine et al. (20071 point out, care should be taken 



when computing dust masses using measured dust temper- 
atures since the thermal emission per unit dust mass at A 
is oc(e'"^/^^ - 1)"^ At wavelengths A <360Atm, dust tem- 
perature has a profound effect on dust mass, since Bi^(r) is 
very dependent on T (so that a 4K difference between 18 
and 22 K results in a 150% increase in Md). At A >450/im, 
Md is less sensitive to dust temperature. Note that the dust 
temperature used in Eq [S] is the T from Eq [3] and not the 
"peak wavelength" dust temperature Td characterizing the 
system as a whole. 

Due to poor constraints on dust absorption coefficients 
and no data points longward of ~450 /im, dust masses for 
the GOALS sample were not reported until the work of [U| 
|et al.| ( [20T2l ), who use the method described in this paper, 
Eq [8j to calculate dust masses at 850/xm using a dust ab- 
sorption coefficient of K850 = 0.15 m^ kg~^ ( [Weingartner fc 
|Draine||2001| [Dunne et al.| |2003). Figure [t] plots the galax- 
ies' dust masses against IR luminosity, which highlights two 
important facts. First, the calculation of dust mass is very 
sensitive to dust temperature, Md oc LirT"^. This high- 
lights that dust mass should only be estimated when dust 
temperature is adequately constrained. The second notice- 
able detail of Figure [7] is the lack of correlation between the 
two physical properties, Lir and Aid- The lack of correla- 
tion is a testament to how large amounts of dust can be 
formed quite quickly in starburst galaxies (in contrast, star- 
forming galaxies with much more moderate star formation 
rates overall have dust masses several orders of magnitude 
smaller, <1O^~^M0). This also points out that scaling rela- 
tions where star formation rate is inferred directly from a 



long wavelength flux density measurement, e.g. ^gso, with- 
out good constraints on Lir or T, are highly uncertain. 



4.4 Future Improvements 

The next decade will see great improvements in our under- 
standing of dust emission in galaxies thanks to vast repos- 
itories of data from large field-of-view mapping bolometers 
and efficient follow-up with interferometric imaging. Facil- 
ities less limited by confusion noise (e.g. ScuBA-2 450/im 
mapping on JCMT, and future 200-850/,im data from the 
Cornell-Caltech Atacama Telescope, CCAT) will make in- 
frared photometric measurements more accurate, precise 
and more numerous. Submillimeter spectrometers (e.g. like 
Z-Spec operated at the Caltech Submillimeter Observa- 
tory and Atacama Pathfinder Experiment, e.g. see |Brad-| 
[ford et al.|[2009| also those being designed for CCAT) will 
measure spectra for individual sources in the infrared and 
submillimeter, providing crucial constraints on emissivity 
while also confirming redshifts through CO lines and assess- 
ing their contribution to continuum photometry (e.g. see 
Smafl et al.|2011[ ). 



With this superb data, infrared SED fits can be signif- 
icantly refined. Better photometric constraints mean that 
more free parameters can be reintroduced into SED fitting, 
providing a more concrete physical context. What does the 
value of P imply for the distribution of cold dust? What 
does the value of a imply for the radial dust distribution? Is 
there a wavelength regime where the assumption of a pow- 
erlaw distribution of warm dust temperatures does not hold, 
and is that in turn related to the dust distribution or pres- 
ence/lack of a luminous AGN? Are any galaxies better fit 
with two independent greybodies rather than one, and does 
this indicate two distinct dust reservoirs not yet mixed? 

Furthermore, interferometers like ALMA will make it 
routinely possible to resolve dust emission on kpc scales 
within distant galaxies (as the Submillimeter Array, SMA, 
has done for some local ULIRGs; [Wilson et al.|[2008| . Spa- 
tial mapping of the dust distribution in multiple infrared 
channels can then be used to assess the mid-IR powerlaw 
inferred dust distribution. Combined with resolved molecu- 
lar gas maps obtained through CO emission, routine tests 



of the Schmidt-Kennicutt relation (Schmidt 1959 Kenni- 
[cutt 1998) can be performed on a variety of gala:xy types at 
different stages of their star formation histories. 

A further refinement which can be made to SED fitting 
described in this paper is the introduction of PAH emission 
and Si absorption modeling on top of the suggested dust 
continuum. Given the wide variety of mid-infrared proper- 
ties for infrared starbursts, this modeling is best fit inde- 
pendently from the mid-infrared powerlaw dust continuum. 
This decoupling of the two will be possible with deep near- 
and mid-IR data from the James Webb Space Telescope 
(JWST). 



5 DISCUSSION & CONCLUSIONS 

The characterization of infrared-bright galaxies has become 
increasingly important in recent years with the introduc- 
tion of new mid- to far-infrared observing facilities, includ- 
ing Herschel Space Observatory and the Atacama Large Mil- 
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limeter Array (ALMA). Infrared SED fitting techniques vary 
widely in the literature, from simple greybody fits to detailed 
dust emission modeling templates. 

Section |2.f | has presented an SED fitting technique 
which can be used to fit a wide range of IR data, from sources 
which have only 3 IR photometric points to sources with 
>10 photometric points. These SED fits do not account for 
PAH emission in the mid-infrared, although they produce 
accurate estimates to a source's integrated IR luminosity, 
dust temperature and dust mass. The fitting technique is 
based on a single dust temperature greybody fit linked to 
a mid-infrared powerlaw, fit simultaneously to data across 
~5-2000/im. Without inclusion of the mid-infrared power- 
law, dust temperatures are overestimated and the short- 
wavelength data (<50/im) fit is likely very poor. From the 
measurement of the mean mid-infrared powerlaw slope a, 
the dust radial density profile in most local (U)LIRGs is 



fairly shallow. 



-1/2 



. IDL code for the fitting procedure 



has been made publicly available. 

This SED fitting procedure is contrasted with fits to in- 
frared template libraries generated through dust-grain mod- 



eling, inclu ding those ofjChary fc Elbaz ( 2001 1, Dale & Helou 
( 20021) and ISiebenmorgen & Kriigell (120071). SED fit quality 



is categorized by goodness of fit to the GOALS local LIRG 



and ULIRG sample photometry reported in U et al. (2012) 



Since the GOALS sources are the most extensively surveyed 
infrared galaxies to date, they provide a good testbed for 
SED fitting reliability. The SED templates and original for- 
mulation of luminosity for the RBGS sample do not fit the 
data as well as the fitting method described in section [27lj 
leading to overestimated IR luminosities by 10%. Similarly, 
dust temperatures estimated for the original RBGS sam- 
ple ( Sanders et al.|2003 l were overestimated by ~10K, due 
to the different fitting method used (single greybody versus 
single greybody plus mid-IR powerlaw). Dust temperatures 
cannot be well constrained with SED templates because they 
are quantized and only peak at certain wavelengths corre- 
sponding to fixed dust temperatures. Clarification is offered 
on the calculation of dust mass, its dependence on dust tem- 
perature, and the lack of correlation between dust mass and 
IR luminosity. 

This comparative study should be useful to highlight 
some of the current problems facing FIR SED fitting tech- 
niques, and the details often overlooked in ambitious anal- 
yses of infrared-galaxy population trends for galaxies both 
near and far. 
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